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Abstract 

We perform lattice Monte Carlo simulations for up to 66 unitary fermions in a finite box us- 
ing a highly improved lattice action for nonrelativistic spin 1/2 fermions. We obtain a value of 
0.3661qq]^^ for the Bertsch parameter, defined as the energy of the unitary Fermi gas measured 
in units of the free gas energy in the thermodynamic limit. In addition, for up to four unitary 
fermions, we compute the spectrum of the lattice theory by exact diagonalization of the transfer 
matrix projected onto irreducible representations of the octahedral group for small to moderate size 
lattices, providing an independent check of our few-body simulation results. We compare our exact 
numerical and simulation results for the spectrum to benchmark studies of other research groups, 
as well as perform an extended analysis of our lattice action improvement scheme, including an 
analysis of the errors associated with higher partial waves and finite temporal discretization. 
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I. INTRODUCTION 

Unitary fermions have gained widespread attention from theorists, particularly since their 
successful creation in experiments involving trapped, ultracold atoms. The universal nature 
of this system promises applications to many fields; for example, it has been suggested as 
an expansion point for an effective field theory for nuclear physics. Numerical studies have 
been essential to the progress in our knowledge about unitary fermions due to the strongly 
coupled nature of the interaction, which renders standard perturbative techniques unreliable. 

Unitarity corresponds to the idealized limit in which the s-wave scattering length becomes 
infinite and the interaction range vanishes, or equivalently, the two-particle s-wave scattering 
phase shift 60 = 7i/2. For a homogeneous system of two-component fermions at unitarity, 
the only relevant scale is the density p = N/V, with A^ = N-^ + N^^. Consequently, the 
ground state energy is related to that of noninteracting fermions by 

E{p)=^EFree{p), (l) 

where Epreeip) = 3NEp{p)/5, Ep{p) = kl/{2M) is the Fermi energy, and kp = {3n^ pY^^ 
the Fermi momentum. The dimensionless parameter ^, known as the Bertsch parameter [T] , 
is of particular interest because it is the unique parameter which relates zero temperature 
thermodynamic quantities between the unitary and free Fermi gas. Several experimental 
groups have measured the Bertsch parameter using a variety techniques involving ultra-cold 
trapped atoms pHTS] . On the theoretical side, in addition to analytical calculations [Ti] - [27] , 
there has also been a substantial number of numerical studies of unitary fermions from the 
microscopic theory using quantum Monte Carlo and other techniques [28] - H^ . Many of these 
studies use the Bertsch parameter as a benchmark calculation. 

At very low energies and large scattering length the detailed structure of the inter-particle 
potential become irrelevant, and the system is well-described by an effective theory for spin 
1/2 fermions with a zero-range contact interaction: 

Here, the field ip = {4'tyi'i) i^ a two-component spinor and the coupling Co is tuned to 
an 0{1) critical value, determined nonperturbatively by an exact, analytic evaluation of 
the two-particle scattering amplitude A^^{p) = p cot 6{p) — ip at zero external momentum 



p= ^/ME |15 



Perhaps the simplest lattice construction of Eq. |2] at finite chemical potential was devel- 
oped by Ref. |17], and employs a nonpropagating real scalar field (p of mass w? = I/Cq to 
induce two-body interactions between fermions of opposite spin through the type of interac- 
tion (pip'^Tp. ^ This lattice construction was shown to be free of the fermion "sign problem" 
at finite density, a problem that is notorious for rendering numerical simulations of certain 
fermionic theories at finite density impractical, due to the presence of a complex effective 
action for the scalar field (f). 

We have since developed a highly improved lattice theory based upon the construction 
of Ref. |17j, allowing us to study up to 70 unitary fermions confined to a harmonic trap 
[50] and up to 66 unitary fermions confined to a finite box. Results of the latter study 
are described in detail in the proceeding sections. Here, we summarize some of the salient 
features of our construction: 

1. We employ open boundary conditions in the time direction, preventing fermion propa- 
gation from wrapping "around the world." This choice eliminates 0-dependence in the 
fermion-determinant obtained upon integrating out the fermion degrees of freedom in 
the path-integral, and therefore yields a trivial effective action for the scalar auxiliary 
field and eliminates the need for importance sampling in the simulation. 

2. Due to our choice of temporal boundary conditions, simulations must be performed at 
zero temperature and zero chemical potential. The energy of systems at finite density 
are obtained by studying the long-time exponential fall-off of multi-fermion correlation 
functions. 

3. We use a continuum single particle dispersion relation for fermions, thus reducing 
lattice discretization errors. 

4. We introduce Galilean invariant derivative interactions which allow us to eliminate 
higher order terms in the effective range expansion for pcot 5o(p). We are thus able to 
simulate fermions close to the unitary limit even at small lattice volumes. 

Several recent papers have indicated that lattice Monte Carlo methods can be affected by 
large systematic errors due to a finite filling factor [^THSB] . For a given number of particles, 

^ This technique is often used in lattice simulations involving quartic fermion interactions, and is commonly 
known as the Hubbard-Stratonovich transformation I481I49I, named after its inventors. 



this systematic error corresponds to a dependence of the Bertsch parameter on the number 
of lattice sites. The improvements referred to above are crucial in reducing these errors. In 
this paper, we provide an extensive discussion of the discretization errors which remain after 
improvement, based upon an analysis of the Symanzik action [5H 155] . 

The organization of this paper is as follows: In Sec. [TT| we summarize our highly improved 
lattice construction for numerically simulating untrapped unitary fermions and provide de- 
tails regarding the construction of multi-fermion correlation functions used to extract the 
ground state energy of the system. In addition, the method used to tune two-body couplings 
to the unitary point is briefly reviewed. The details of the lattice construction are discussed 

we 



at greater length in Ref. [50]; here we only provide the main ingredients. In Sec. Ill, 
present exact spectrum results for the two- and three-fermion systems on a lattice at finite 
volume, and use those results to try to understand the systematic errors associated with 
our construction due to interactions from higher partial waves and temporal discretization 
errors; a description of how the multi-fermion transfer matrix is constructed is provided 
in the appendix, along with the construction of projection operators onto the irreducible 



representations (irreps) of the octahedral group. In Sec. IV, we summarize the techniques 
used for extracting the energies of up to 66 unitary fermions in a finite box, and present 
simulation results for the few- and many-body system, including an estimate for the Bertsch 
parameter. 

II. LATTICE CONSTRUCTION 

A. Action 

We consider a highly improved lattice theory for an interacting system of nonrelativistic 
spin 1/2 fermions of mass M on a T x L^ Euclidean space-time lattice with temporal extent 
T and spatial extent L. The sites of the lattice are labeled by integers r G [0, T) in the 
time direction and Xj G [0,L) in the space directions with j = (1,2,3). Throughout this 
work we impose open boundary conditions in the time direction and periodic boundary 
conditions in the space directions. Unless otherwise noted, we measure all quantities with 
dimensions of energy in units of the inverse temporal lattice constant 6^^ and all quantities 
with dimensions of length in units of the spatial lattice constant bs- 



The lattice action for this theory is given by |50j : 

S = Y,i^lKi;^, (3) 

where V'o- and ip'l are Grassmann valued (T x L^)-diniensional vectors associated with each 
spin component a = (f, \.), and i^ is a (T x L^)-diinensional matrix of commuting numbers. 
The matrix elements of the fermion operator K are given by 

where 

f eP'/2^ , IpI < A 
Dp,p' = 5p,p, X < , , (5) 

I oo , IpI > A 

and 

Xp,p.(r) = Vp' + C'^'iP - P')0P-P'(r) . (6) 

The matrix elements are labeled by a time coordinate r and by a three-momentum pj = 
27mj/L, where Uj G [— L/2,L/2) for a periodic spatial lattice (assuming even L). 

Two-body interactions are induced by the periodic field 0p(t), defined as the spatial 
Fourier transform of a random auxiliary field (pxi^) in position space, which satisfies the 
conditions: 

(0x(r)) = , (0x(r)0x'(r')) = K^'Sry ■ (7) 

Throughout this work, we take 0x(t) to be a Z2-valued field with probability distribution 
p(0) = {6i^ff, + 5_i^(^)/2 for all x and r. ^ The two-body coupling C(p) is a periodic function 
of momenta and is given by the operator expansion: 

No-l 

M 



A-jr 

^(P) =m11 C^nOUv) , (8) 



n=0 

up to some fixed order No ~ 1- Throughout this work we use the operator basis: 



C 



2n 



;p) = Mo"(l-e-P'/^^")" , (9) 



^ Although we use Z2 auxihary fields in this work, Gaussian distributed (/) fields with probability distribution 
p{(^") — e^'^ ' ^ for every x and r would work equally well. 



where p is taken to be a periodic function of p and satisfies the relation p^ = p^^(A — 
IpI) + A^^(|p| — A) for p G BZ^ with BZ denoting the first Brillouin zone. Note that at 
low momenta, the operators satisfy the low energy expansion 02n(p) = P^" [1 + C^(p^)]) 
irrespective of the mass parameter Mq. For simplicity, we choose M = Mq, although these 
mass parameters need not be the same. 

The partition function for the lattice theory is obtained by integrating out the fermionic 
degrees of freedom, yielding a path-integral over the auxiliary fields given hj Z = 
J[d(l)]p{(j))detK. Due to the upper block tri-diagonal form of the fermion operator, one 
may show that the fermion determinant is given by detK = detD^, and is independent of 
the auxiliary field. Importance sampling for the field in a Monte Carlo simulation thus 
reduces to generating random field configurations distributed according to the trivial distri- 
bution p(0). Hence, a full simulation of the theory is the same as a quenched simulation. 

Upon integrating out the fermion degrees of freedom, expectation values of operators 
involving the fermion fields ip^ and ip^, such as multi-fermion correlation functions, reduce 
to the expectation values of appropriately contracted fermion propagators K^^. A fermion 
propagator, which evolves a single fermion state from time slice zero to time slice r over a 
given background auxiliary field, may be expressed using the simple recursive formula: 

K-\t; 0) = D-^X{t - 1)K-\t - 1; 0) (10) 

with K-\0;0) = D-\ 

B. Observables 

Multi-fermion correlation functions are constructed using sources formed from a direct 
product of single particle states jaj, o") labeled by spin a = {f, l) and wavefunction quantum 
numbers ctj, where i = 1 . . . ,N. Separable sources of this form are typically favored in 
Monte Carlo simulations due to the nature of the algorithms used. In this study, we choose 
eigenstates of the noninteracting Hamiltonian (i.e., a = p) as our single particle states, and 
use the free iV-fermion ground state as our source by filling states in momentum space up to 
the Fermi surface. Specifically, we consider an initial state of the form |iV/2,'|') ^'P\N/2,D, 
where 

\N,a) = ei,,..,,i^|pii,(T)® ...® |pi^,o-) , (11) 
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and P is a parity flip operator satisfying V\pi,o-) = \—pi,a). The specific momenta pi used 
in this construction are tabulated in Table [Tj Note that the parity flip operator is necessary 
to ensure that the total momentum of the initial state vanishes when N/2 is not equal to 
one of the closed shell values (e.g., 1, 7, 19, 27, 33, ... ). 

Although the form of source wavefunctions is constrained by separability due to the nature 
of our algorithm, greater freedom is allowed for the construction of sink wavefunctions. In 
order to maximize the overlap with the unitary Fermi gas ground state, sink wavefunctions 
are constructed by considering a direct product of N/2 spin-paired two particle states, 
I^P) (S) . . . (8) I^E'), following the approach of Refs. [2H1 ED]- The two-particle states used in 
such a construction are given by 



^^ = vT. ^(p)ip.t)®i-P,;) 



V 

peBZ 

with the two-particle pairing wave-function given by 

„-;3|p| 
^(P) 



(12) 




(13) 



The functional form of the pairing wave function \^(p) in momentum space was obtained 
by considering the Fourier transform of the continuum pairing wavefunction \l/(x,.ei) = 
g-/3 l'''-e'l/|xre;| for two fermions in the unitary regime in position space, where Xrei = x-f- — xj, 
is the relative coordinate between the pair. Note that in the continuum and away from uni- 
tarity, P is identified with the scattering length a. Here, we introduce P and \l'o, the zero 
momentum component of the wave function, as tunable free parameters which may be varied 
in order to improve the overlap with the ground-state wave function; the specific values used 
in our simulations will be discussed in Sec. HVl 

Given the sources and sinks defined above, the multi-fermion correlation function for 
A^''^ = N-^ = N/2 fermions used in our simulations is finally given by: 

C^(r) = (det^^^(r)), (14) 

where 

sfM^) = E ^(q)(q|^"'(r,o)|p„;)(-q|i^-i(r,o)|p„t) , (15) 

q&BZ 

and S"^^ is an A^/2-dimensional matrix with indices i,j = 1, . . . , N/2. Note that the deter- 



minant appearing in Eq. [14] ensures that fermions of the same species are properly antisym- 
metrized. 

C. Parameter tuning 

The lattice action defined in Eq. Is] contains Nq couplings C2n {n = 0, . . . , Nq — 1) which 
must be tuned to scattering data. The details of our tuning procedure are described at 
length in Ref. [50], and is similar to the method used in Ref. [56]; here we summarize the 
main points. The couplings C2n are tuned by matching the lowest No s-wave eigenvalues 
A = e~^ of the two-body transfer matrix defined on the lattice at finite volume onto the 
lowest No solutions to Liischer's formula, given by [571460] : 



1 f'pL\^ 

pcot6o{p) = -—S{7]) , r]= ( — ] , (16) 

TvL \2t^ J 



where S{ri) is the three-dimensional Zeta function: 

v^ 1 



S{'q) = lim 

A— !>oo 



/ T^9 47rA 

|j|<A l-JI ' 



(17) 



In the unitary limit, the solutions to Liischer's formula are just the roots of the function S{ri), 
which may be easily calculated numerically and are related to the energies hj p = y~ME. It 
was shown in Ref. [SO] that this tuning procedure may be used to systematically eliminate 
the leading No terms in the effective range expansion 

pcot 5o(p) = -- + - ^ r„_ip2" (18) 

n=\ 



up to negligible residual contributions to r„_i for n < Nt, 



o- 



III. EXACT RESULTS FOR FEW-BODY STATES 

In Appendix |X] we derive an exact expression for the A^-particle transfer matrix, as well 
as projection operators onto the center of mass (cm.) frame and the irreducible represen- 
tations (irreps) of the octahedral group Oh- Using these results, we have performed exact 
diagonalization of the A^= 1 + 1, A^ = 2 + l and A^ = 2 + 2 unitary fermion transfer ma- 
trices on small to moderate lattice volumes. ^ Armed with exact numerical results for the 
eigenstates and energies, we investigate the systematic errors associated with partial wave 
scattering from nonzero angular momentum interactions in the two-body sector, as well as 
finite volume and lattice spacing artifacts in the three-body sector. In addition to these 
studies, we perform consistency checks with our numerical simulation in both the three- and 



four-body sectors. Those results, however, will be presented in Sec. IV 



A. Two unitary fermions 



As discussed in Sec. II C the unitary limit corresponds to the limit that pcot(5o = 
for all momenta, where (5o is the s-wave scattering phase shift. One tacitly assumes that in 
addition to this limit, the effects of scattering in higher partial waves are negligible compared 
to those of s-wave scattering. The latter is a condition that becomes arbitrarily valid at low 
energies (or equivalently, in the dilute limit), when s-wave scattering becomes the dominant 
contribution to the total scattering cross-section. 



^ The dimensionality of the transfer matrix after projecting onto the c.m. frame and irrep r scales roughly 
like (AL/27r) /48. For a fixed amount of computer memory or computing time, the maximum 

allowable lattice size decreases sharply with increasing N . 



10 



In its simplest form, our lattice action involves only a single four-fermion contact inter- 
action (i.e., taking Nq = 1 in Eq. [s]), which allows for scattering in s- waves, but not in 
the higher partial waves. In this case, there is only a single parameter Cq, which may be 
used to tune the s-wave scattering length to infinity. However, since the effective range and 
higher order shape parameters appearing in the effective range expansion for p cot 5q are 
typically 0{1) in lattice units, extrapolations in the lattice volume are required to elimi- 
nate systematic errors associated with those parameters. One may reduce the systematic 
errors by introducing higher derivative lattice operators (i.e., taking Nq > 1) as described 



in Sec. II C In doing so, however, one in turn introduces systematic errors associated with 
higher partial-wave interactions which are attributed to the fact that the lattice interactions 
no longer involve fermions at a single lattice site, but also all the neighboring sites as well. 

In light of these considerations, it is important for us to estimate the size of systematic 
effects attributed to higher partial wave interactions in our lattice theory when Nq > 1. 
We proceed by studying the energy eigenvalues associated with higher partial-waves in the 
two-body sector, and particularly study their deviation from the energies expected for two 
noninteracting fermions. It is well known that the eigenstates of the lattice theory transform 
as irreps of the octahedral group, and that those irreps contain specific angular momentum 
components in the infinite volume and continuum limits [6lj. The decomposition of orbital 
angular momentum eigenstates into the various irreps of the hypercubic group are provided 
for reference in Table [IT} Note that even i correspond to the positive parity irreps, whereas 
odd i correspond to negative parity irreps. 

A cursory examination of Table O shows that the Af lattice eigenstates contain angular 
momenta components £ = 0, 4, 6, 8, ... in the continuum and infinite volume limits. Similarly, 
the Tf eigenstates possess the components i = 1,5,7,9,.... Generally speaking, the size 
of the effects of interactions with angular momenta £ < 6 and i = 9 may be deduced by 
studying the energies of the two-body lattice eigenstates classified by their irreps under Oh. 
By studying the deviations in the lattice eigenvalues t], defined by —log A = jj (^) i], 
where A is an eigenvalue of the two-body transfer matrix, in a given irrep from those of 
noninteracting fermions r/*, which take integer values, we may estimate the size of the 
effects from scattering in higher partial waves. Figures [I] and |2] show the deviation r^/?]* — 1 
as a function of rj for the entire spectrum for unitary fermions on an L = 8 lattice with 
1 < Nq < 4 and an L = 16 lattice with 2 < A^^ < 5. We find that the largest deviation is 

11 



TABLE II: Decomposition of angular momentum eigenstates into irreps of Oh- Also indicated are 
the Oh irreps containing i as their lowest-lying state (LLS). 
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approximately 2% in the smallest box size considered, while the deviations are considerably 
smaller for the L = 16 box (< 0.6%). 

Using the generalization of Liischer's formula for s-wave scattering, we may determine 
the scattering phase shifts for the higher partial waves. For p-wave scattering, if one assumes 
tan54 ^ tan^i, one finds (see, for example, [62]): 



p'^ cot 61 (p) 



27r 

T 



^2VS{V) 



(19) 



where rj and S{ri) are defined in Eq. 16 Plugging the lattice eigenvalues obtained for the 



Tj" irrep shown in Fig. ^into the right-hand side of Eq. 19, we obtain a lattice prediction 
for p'^cot^i. In Fig. [3| the scattering phase shift 5i obtained by this procedure is plotted 
as a function of rj for L = 8 and L = 16, and for Nq 



1,2,3,4 and 5. For reference, also 



shown in this figure is the scattering phase shift 6q obtained from Eq. 16 



More generally, one can show from the results of Ref. [62] that if 7] is sufficiently close to 



T]*, then for the partial waves 



Siip) 



0, 1, 2, 3, 4, 5, 6 and 9 one finds: 



iv' 



x3/2 



9e{v*] 



iri/ri*-l) + 0{ri/ri*-lf 



(20) 



where geirj*) is some non-zero calculable numerical factor. For s- and p-waves, goif]*) = 
giijf) = (i(?7*)/(27r^), where d{ri*) equals the number of times the integer triplet j satisfies 
IJI = Tj* for a given pole rj* (i.e., an integer taking the value: 1, 6, 8, 12, 24, or 48 depending on 
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FIG. 1: Energy spectrum, given by —log A = jj {~^) V^ of ^^o unitary fermions of mass M = 5 
in a finite box of size L = 8 for up to four tuned couplings; rf are the corresponding eigenvalues 
of the noninteracting theory. Eigenvalues are labeled by dimensionality of the irrep to which they 
belong: circle (A), square (E), and triangle (T), as well as color coded according to the lowest 
orbital angular momentum component contained in each irrep: red {I = 1), orange {i = 2), yellow 
(£ = 3), green (^ = 4), cyan (^ = 5), blue (^ = 6), and violet (f = 9). 

7]*)] for £ > 1, the expression for ge{ri*) is more complicated (involving spherical harmonics), 
and therefore is not provided here. From Fig. [T] and Fig. [2| it is evident that the deviation 



of higher angular momentum modes diminish with i, and likewise based on Eq. 20 so must 
the corresponding phase shifts. 



B. Three unitary fermions 

Here we present results from exact numerical diagonalization of the A^ = 2 + 1 unitary 
fermion transfer matrix in an effort to better understand the effectiveness of our parameter 
tuning method. In addition, we use the exact numerical results to investigate the finite 
spatial and temporal discretization errors of the lattice theory. Exact numerical results for 
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FIG. 2: Energy spectrum, given by —log A = -j-^ {~^) V^ of ^^o unitary fermions of mass M = 5 
in a finite box of size L = 16 for up to five tuned couplings; r/* are the corresponding eigenvalues 
of the noninteracting theory. Eigenvalues labeled according to irrep, as described in Fig. [T] In the 
case iVo = 1 (shown in Fig. [llfor L = 8, but omitted here), the deviation 77/77* — 1 is exactly zero 
for all nontrivial irreps. 

the A^ = 2 + 1 system were obtained for the A^ irrep for L = 4, 6, 8 and 10 and for up to 
Nq = 5 tuned couphngs. At small lattice volumes, we were limited to fewer than five tuned 
couplings for reasons discussed in Ref. [50j. 

Fig.Hshows the energies of the ground and first excited states in the Af irrep as a function 
of volume, along with the extrapolated continuum limit, infinite- volume results obtained in 
Ref. [63]. Note that these energies are measured in units of the non- interacting energy for 
three zero-momentum fermions, Epree = (2,71 / L)'^ / M . Extrapolation of the No = 1 results 
using a linear fit in L^^ to the data yields better than 1% agreement with results from 
Ref. ^5; deviations may be attributed to the long extrapolation performed on our exact 
finite volume energies. Note that for A^^ > 1 tuned couplings, the lattice energies show 
substantially improved agreement with exact results for the ground state energy-even at 
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FIG. 3: 6o and 6i as a function of the diniensionless parameter rj for unitary fermions of mass 
M = 5 on an L = 8 and L = 16 lattice, obtained from Eq. 16 and Eq. 19, respectively. Dashed 
lines correspond to the phases 7r/2 (left) and zero (right). 
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FIG. 4: Left: ground-state energy of A^ = 2 + 1 unitary fermions of mass M = 5 in the A^ irrep 
as a function of 1/L for up to four tuned couplings. Right: first excited state energy of A^ = 2 + 1 
unitary fermions of mass M = 5 in the A^ irrep as a function of 1/L for up to four tuned couplings. 
Dashed lines in both plots correspond to the infinite volume exact results of Ref. [63]; solid lines 
are infinite volume extrapolations based on a linear fit to all Nq = 1 data (deviations in the 
extrapolation reflect systematic effects of unaccounted for higher order corrections). 
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FIG. 5: Ground-state energy of A^ = 2 + 1 unitary fermions in the Af irrep as a function of 1/M 
for one (left) and four (right) tuned coupHngs with L = 8. Blue dashed line corresponds to a fourth 
order polynomial fit in 1/M; black dashed line indicates the mass value at which all simulations 
have been performed. 

small to moderate volumes-compared to the No = 1 results. The improvement with Nq 
is less pronounced in the excited Af state, however, for which we do not have a rigorous 
explanation. 

Although tuning more than one s-wave two-body operator results in increased improve- 
ment in ground state energies for L = 4 and 6, little improvement in energies is evident 
among No = 2,3,4 and 5 results at L = 8 and L = 10. It is possible that this peculiar 
behavior may be due to the effects of an untuned £ = 1 two-body operator (giving rise to 
L~^ scaling) or three-body operators (which are expected to contribute starting at L~^'^^). 
On the other hand, the leading volume corrections from untuned two-body i = operators 
scale as l/L^^"^"^. At large L, the former volume corrections may dominate irrespective of 
No > 1, whereas for small L the latter corrections becomes non- negligible even for large 
No, giving rise to stronger A^o-dependence of the energies at small L. 

In Fig. |5} we show the fermion mass-dependence of the A^ = 2 + 1 ground state energy for 
A'c' = 1 aiid No = 4 on an L = 8 lattice. The data was fit using a fourth-order polynomial 
in 1/M and extrapolated to the M — )■ oo limit. Since the physical (dimensionful) mass 
is equal to Mbr/b^, taking M — )■ oo is equivalent to taking the temporal continuum limit 
6t- — ^ while keeping the spatial lattice spacing bs and physical mass held constant (in other 
words, the lattice mass parameter may be viewed as an anisotropy factor). From Fig. [51 we 
see that the temporal discretization effects on the three fermions system at M = 5, the mass 
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TABLE III: Simulation parameters for untrapped fermions (M = 5) using the pairing wave function 

100. For ensembles denoted by an asterisk, we have used ^{p) = 

^conf is the total size of the 
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with ^ 







given by Eq. 

l/(/3 — e~P /(2A2)) for the pairing wave function rather than Eq. 
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ensemble 


observable data was 


block-averaged in blocks of size Nconf/^B prior 


to analysis. 


N 


L 


T 


No 




Ncor,f 


Nb 


4 


4 


24 




0.1 / 0.4 


190M / 330M 


190 / 330 


4 


6 


36 




0.2 / 0.35 / 0.5 


300M / 300M / 300M 


300 / 300 / 300 


4 


8 


48 




0.15 / 0.25 


500M'' / 500M 


250 / 250 


4 


10 


64 




0.25 / 0.3 / 0.32 


400M / 400M / 400M 


400 / 400 / 400 


4 


12 


54 




0.25 / 0.31 / 0.35 


400M / 400M / 400M 


400 / 400 / 400 


4 


14 


54 




0.37 / 0.43 


680M / 680M 


340 / 340 


4 


4 


24 


2 


1.01* / 0.1 


2.67B / 1.29B 


300 / 300 


4 


6 


36 


2 


0.1 


200M 


200 


4 


10 


42 


2 


0.15 / 0.19 


400M / 300M 


400 / 300 


4 


8 


48 


5 


0.05 / 0.15 


llOM / llOM 


110 / 110 


4 


10 


64 


5 


1.01* / 0.07 / 0.1 


150M / lOOM / lOOM 


300 / 200 /200 


4 


12 


64 


5 


0.15 / 0.2 / 0.25 


150M / 60M / 180M 


300 / 120 / 360 


4 


14 


64 


5 


0.2 / 0.25 / 0.3 


250M / 250M / 130M 


250 / 250 / 130 


4 


16 


64 


5 


0.25 / 0.3 / 0.35 


390M / 390M / 390M 


390 / 390 / 390 


4 


18 


64 


5 


0.25 / 0.3 


350M / 350M 


350 / 350 


< 66 


10 


54 


5 


1.0 


40M 


400 


< 66 


12 


64 


5 


0.5 


40M 


400 


< 66 


12 


54 


5 


0.75 


40M 


400 


< 66 


12 


64 


5 


1.0 


40M 


400 


< 66 


14 


36 


5 


0.5 


19M 


190 


< 66 


14 


36 


5 


0.75 


35M 


350 


< 66 


14 


72 


5 


1.0 


39M 


390 


< 66 


16 


54 


5 


0.9 


20M 


200 


< 66 


16 


30 


5 


0.9 


20M 


200 


< 66 


16 


30 


5 


0.75 


63M 


315 



"For this ensemble, we have used \E'o — 50. 

value used in our few- and many-body simulations, is roughly 0.5% for A^^ = 1 and 0.1% for 
No = 4 tuned operators. Time discretization errors are therefore likely negligible compared 
to other systematic and statistical uncertainties in our few- and many-body simulations. 
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TABLE IV: Simulation parameters (i.e., tuned couplings) for untrapped unitary fermions (M = 5) 
for various lattice volumes and Nq values. 



L 


Co 


C2 


Q 


Ce 


Cg 


4 


0.673068 


- 


- 


- 


- 


6 


0.689184 


- 


- 


- 


- 


8 


0.680971 


- 


- 


- 


- 


10 


0.684858 


- 


- 


- 


- 


12 


0.679787 


- 


- 


- 


- 


14 


0.684345 


- 


- 


- 


- 


4 


0.333477 


0.1552055 


- 


- 


- 


6 


0.428091 


0.1128065 


- 


- 


- 


10 


0.455289 


0.0939424 


- 


- 


- 


8 


0.931735 


-2.1243485 


2.2200002 


-0.7798253 


0.08856646 


10 


0.585273 


-0.1507720 


0.2120923 


0.0974153 


0.01455297 


12 


0.544064 


-0.0354881 


0.0770458 


-0.0433194 


0.00783886 


14 


0.547526 


-0.0218146 


0.0489023 


-0.0291435 


0.00588451 


16 


0.537953 


-0.0042753 


0.0284083 


-0.0211698 


0.00492156 


18 


0.547918 


-0.0111534 


0.0375212 


-0.0279792 


0.00613497 



IV. SIMULATION RESULTS FOR FEW- AND MANY-BODY STATES 

A. Ensembles and parameters 

Our numerical studies of untrapped unitary fermions consisted of two parts: 1) high 
precision calculations for A^ < 4 fermions, intended for investigating the systematic errors 
associated with finite volume artifacts in the few-body system, and 2) simulations for up to 
A^ = 66 fermions in order to extract a thermodynamic limit value for the Bertsch parameter. 



Ensemble details for each of these studies are provided in Table III the C2„ values used for 



a given L and Nq at M = 5 are provided in Table IV Throughout our studies, we have used 



an ultra-violet cutoff of A = 0.999997r. For our few-body studies we generated ensembles 
of size Nconf ~ lOOM-lB on lattices ranging from L = 4 — 18 in size and No = 1, 2, and 
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5 tuned couplings. The total CPU time required for all fifteen few-body simulations was 
approximately lOOK CPU-hours. In our many-body studies of up to 66 unitary fermions, 
we generated ensembles of size N^onf ~ 20M-60M on lattices ranging from L = 10 — 16 and 
with No = 5 tuned couplings. A total of approximately 450K CPU-hours was required to 
generate all ten many-body ensembles. 

Observables were measured on each ensemble and the results were averaged into Nq blocks 
of size Nconf/Nis prior to analysis. For even A^, correlation functions were measured using 
Eq. [14| and appeared to be insensitive to the free parameter \l/o in the sink wave-function 
provided that the parameter is sufficiently large. Throughout this work we therefore fixed 



\I^o = 100 in Eq. [I3]but considered multiple values for the free parameter /3 on each ensemble. 
Details regarding the ensembles and correlation functions used for the case when A^ = 3 are 
provided in j50] 

The lattice action possesses one additional free parameter, the mass term M, which 
controls the anisotropy of the lattice (i.e., a conversion factor between space and time). 
The temporal discretization errors in the many-body problem are controlled by the quantity 
kp/M ~ (N^^^/Ly/M. For fixed br, the temporal discretization errors are subleading in 
the density compared to spatial discretization errors, which are controlled by kp ~ N^^^/L. 
The temporal discretizations are therefore expected to be under control provided M is larger 
than 0{1) in lattice units. Since the decay rate of correlation functions are proportional to 
1/M, we may obtain earlier plateaus in r by decreasing M, however, this comes at the cost 
of increased temporal discretization errors. With these considerations in mind, we find that 
M = 5 is an ideal compromise and use this value throughout all of our studies. 

B. Analysis technique 

Throughout this work, we study the behavior of multifermion correlation functions at 
late time in order to extract information about the low-lying spectrum of the system at 
unitarity. Specifically, correlators have the late-time behavior: 

C(r) = Zoe~'^°'^ -|- excited state contributions , (21) 

where Zo is a complex number which quantifies the overlap between our source and sink 
wave functions and the multi-fermion ground state, and Eq is the corresponding ground state 
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energy of the system. For small numbers of fermions, we may use a conventional approach 
for extracting Eq by studying the plateau region of the effective mass, defined as: 

meffir) = —- log ^, , , . (22) 

The correlator C(r) is typically estimated from an ensemble average of correlators measured 
on random background field configurations and At is an arbitrary positive integer, usually 
chosen to be unity. In the late time limit one finds me//(T) ^ Eq up to corrections that 
are exponentially small in the energy splittings. Using more sophisticated analysis methods, 
one may extract excited states from the effective mass as well. 

For larger A^, however, correlators measured using our simulation algorithm generally 
possess a distribution overlap problem, rendering conventional estimates of the effective 
mass unreliable. The problem is particularly severe when the number of configurations 
is less than on the order of e^^^^^^''^'^''^'^'^'^\ where Epree{N) is the free gas energy of A^ 
fermions [5^. For small numbers of fermions this problem may be overcome with brute 
force by generating very large ensembles, but for large numbers of fermions, brute force 
becomes impractical and one must resort to alternative techniques for reliably estimating 
the effective mass. 

The approach we take for estimating effective masses for large numbers of fermions 
exploits properties of the distribution for multifermion correlators. Particularly, we have 
demonstrated through numerical studies as well as a mean-field calculation that the correla- 
tor distribution function is log-normal in character, thus motivating a method for extracting 
effective masses based on the properties of cumulant expansions. Defined in terms of a 
cumulant expansion, the effective mass is given by 



^^ = ^11^ ['^"(^) - '^"(^ + ^^)] ' (23) 



n=l 



where K„(r) in the nth cumulant of the distribution of the logarithm of the correlator. For 
perfectly log-normal distributed correlation functions, the above expansion truncates exactly 
at second order, whereas for distributions that deviate from log-normal, such deviations are 
encoded in nonzero but presumably small contributions to the sum at order n > 2. In 
practice, the cumulants are estimated from the moments of the logarithm of the correlation 
function, and one must carefully study the effective masses as a function of the truncation 
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order A^^ to determine the ideal value at which statistical errors in the estimate of fi;„(r) are 
comparable to the systematic error associated with the truncation of the expansion. 

In each of our our studies, the ground state energy of the multifermion system has been 
estimated by performing correlated x^ fits of the effective masses to a constant over the 
plateau region at late Euclidean time. In most cases, we considered either two or three 
values of the sink parameter /3 for each ensemble in order to gauge possible systematic 
errors associated with excited state contamination, which may arise due to poor overlap 
with the ground state for a given interpolating operator or due to energy splittings which 
are smaller than the typical inverse time considered. Simultaneous correlated x^ fits to the 
correlation functions were performed using all available correlators on a given ensemble. 

In cases where a plateau in the effective mass plot failed to appear before the onset of 
noise, we instead fit the ground plus excited state using a constant plus an exponential fit 
function. In all cases, the statistical uncertainties were obtained by resampling data using 
the bootstrapping method. In order to take into account systematic effects due to temporal 
correlations and excited state contamination, we varied the end points of the fitting region 
by up to ±3 time steps and regarded the maximum and minimum fit values as our fitting 
systematic errors. The total fitting uncertainty is determined by combining both statistical 
and systematic errors in quadrature. 

C. Few-body Results 

Numerical simulations of A^ = 2+1 and A^ = 2+2 untrapped unitary fermions at zero total 
momentum were performed in order to study finite volume effects as a function of 1/L, as 
well as to make direct comparisons with precision benchmark results of previously reported 
studies. In addition, a comparison with exact diagonalization results of the A^ = 2 + 1 and 
A^ = 2 + 2 unitary fermion transfer matrices on small volumes provide a nontrivial check for 



our lattice simulations. As was the case in Sec. |IIIB[ all few-body energies are measured in 
units of the non-interacting few-body energies, e.g., Ep^ee = (27r/L)^/M for both A^ = 2 + 1 
and A^ = 2 + 2 fermions at zero total momentum. 

In Fig. |6| we plot simulation results for the energy of three unitary fermions at zero total 
momentum on lattice sizes up to L = 16 and for No = 5. These results were originally 
reported in Ref. [50]. Exact ground state energies for the Af irrep obtained from diagonal- 
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FIG. 6: Energy of A^ = 2 + 1 unitary fermions in a zero total momentum eigenstate as a function 
of 1/L^. Blue data points and associated error bars were obtained from numerical simulation, 
short blue dashed lines at L = 8 and L = 10 indicate results from exact diagonalization of the 
three fermion transfer matrix. Red error band indicates the infinite volume extrapolation result 
previously reported in Ref. |50] using simulation data. Black dashed line indicated the exact 
infinite volume result of Pricoupenko and Castin reported in Ref. |63j . 



izing the transfer matrix at L = 8 and L = 10 are indicated in the figure, and agree with 
our simulation results to within errors of 0.16% and 0.18%, respectively. As discussed in 
Ref. |SD] , the leading volume-dependent corrections to the energy of more than two unitary 
fermions with many two-body s-wave operators tuned is expected to be of order L~^, com- 
ing from an untuned two-derivative two-body p-wave operator. Subleading corrections are 
expected to be of order L~^'^^, due to the lowest dimension three-body operator, which has 
i = and scaling dimension 4.67 [USHUH] . Performing a fit to the data using the functional 



form Co -|- Ci/L^ yields an infinite volume extrapolation result of E/Epree = 0-3735lQ;QQQy, 
and is consistent with the exact infinite volume result of Pricoupenko and Castin [63] within 
0.3% uncertainties. 

In Fig. [7], we have summarized simulation results for the ground state energy of four 
unitary fermions for up to A'^^ = 5 tuned couplings and lattice sizes up to L = 18. Exact 
lattice energies obtained for L = 4 are plotted in Fig. [7] for A^^ = 1 and 2 couplings tuned 
to unitarity. In each case the exact ground state energies obtained from the transfer matrix 
are consistent with the simulation results within uncertainties. In a high precision check, 
we found that the ground state energy of four unitary fermions at L = 4 and No = 2 
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FIG. 7: Ground state energy of A^ = 2 + 2 unitary fermions as a function of 1/L. Error bars 
include statistical and fitting systematic errors combined in quadrature. The blue and yellow 
bands represent fit results to Nq = 1 and Nq = 5 data as discussed in the text, with error 
bands reflecting both statistical and systematic errors. Black dashed lines indicate the error band 
obtained from an infinite volume extrapolation of exact benchmark calculations reported in Ref. 
[69] . Short dashed lines at L = 4 indicate energies obtained by exact diagonalizing the four-body 
transfer matrix. 



obtained from ensembles of approximately 45 configurations agreed with exact results to 
within errors of 0.05%. 

For No = 1, the leading volume correction to the ground state energy for four fermions 
will be of order 1/L, due to the untuned effective range operator. To extract the ground 
state energy at L = 00, we therefore used Cq + Ci/L as our fit function for the extrapolation. 
We take into account systematic errors in the infinite volume extrapolation by varying the 
fit interval from L = [4, 14] to L = [10, 14], and obtain E/Epree = 0.2122(40) for the ground 
state energy. For the highly tuned No = 5 case, we expect the leading volume dependence 
for four fermions to be L~^, using the same reasons as for three unitary fermions. Unlike the 
case for three fermions, however, the lowest dimension three-fermion operator is expected to 
have i = 1 and scaling dimension 4.27 rather than i = 0. The reason is that three of the four 
fermions are not restricted to a specific angular momentum state. The subleading volume 
dependence is therefore expected to scale as L~^'^^. Additional subleading terms scale as 
L~^'^^ and L~^ corresponding to the i = three-body operator and the four-derivative p- 
wave and d-wave two-body operators, respectively. By considering the leading L-dependence 
induced by these operators, we use the fit function: cq + Ci/L^ + C2/ Li^'^^ to extrapolate the 
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energy in the case of No = 5. The fit result over the interval L = [10, 16] is shown in 
Fig. [tI and at infinite volume we obtain E/Epree = 0.2130(26). Both our A^^ = 1 and 
Nq = 5 results for the ground state energy of four unitary fermions are consistent with the 
benchmark calculation reported ^69j, within the given uncertainties. 

D. Many-body Results 

To determine the Bertsch parameter, we calculate the ground state energies of up to 
66 untrapped and unpolarized unitary fermions using the many-body ensembles described 



in Table |III[ For small A^, the distribution overlap problem is absent and the conventional 
effective mass defined in Eq. |22] typically shows an acceptable plateau. On the other hand, 
for large A^ the conventional effective mass exhibits a significant overlap problem and we 
generally fail to find plateaus. Examples of each of these scenarios are shown in Fig. [8] 
(upper-left) and Fig. |9] (upper-left) for A^ = 10 and A^ = 50, respectively. The conventional 
effective mass for A^ = 10 in Fig. |8] shows a plateau beginning at around r ~ 23 and the 
ground state energy may be calculated by performing a constant fit to the plateau region 
before the onset of severe noise at r ~ 37. However, the conventional effective mass for 
A^ = 50 in Fig. [9] drifts upward beginning at r ~ 7 and exhibits no plateau before the onset 
of an overlap problem. 



As discussed in Sec. |IVB[ we have used a cumulant expansion technique to overcome 
the distribution overlap problem for large A^. In particular, using the effective mass defined 
in Eq. 23 at moderate truncation orders A^^, we find that rn^r'i {t) exhibits clean plateaus 



for A^ < 66. By performing a constant fit to the plateau region, we estimate the energy 
for each N^^. For small A^, we may verify the cumulant method by comparing energies 
with those measured using the conventional analysis. In Fig. 8, we plot 'm),r'l for A^ = 10 
with Ni^ = 3,5,7 along with the conventional effective mass for comparison. The cumulant 
effective mass with A^^ = 3 shows a clean signal, but does not exhibit a plateau for any given 
time extent. However, as we increase the truncation order A^^, fneft appears to converge, 
and beyond A^^ = 6 the addition of higher cumulants only increases statistical noise without 
changing the plateau. A plot of the energies obtained at each truncation order is presented 



in Fig. 10 (left), along with the energy measured from the conventional effective mass. Our 



best estimation of the energy using the cumulant expansion is obtained for A^^ = 6 using 
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FIG. 8: Conventional (?n,e//) and cumulant {rn^ff) effective mass plots for A^ = 10 unitary 
fermions with At = 2 on a L = 12 lattice. Upper-left panel shows conventional, upper-right shows 
cumulant with N,^ = 3, lower-left shows cumulant with N,^ = 5, and lower-right shows cumulant 
with N,^ = 7. The purple band in the effective mass plots represent fits results to the plateau 
region when one exists; the gray data in the cumulant effective mass plots represent the effective 
mass obtained by using the conventional method. 



the convergence criteria outlined in Ref. [M], and is consistent with the energy obtained by 
the conventional approach. 

For large numbers of unitary fermions, the conventional method fails to exhibit a plateau 
due to the onset of an overlap problem, as demonstrated in Fig. |9] for A^ = 50. We must 
therefore rely entirely on the cumulant expansion to estimate energies in this case. In Fig. [9] 
we plot the effective masses from the cumulant expansion at truncation orders A^^ = 2, 3, 
and 4. We find that irif^r't for small A^^ has a clean signal and the fit results to the plateau 



region shown in Fig. 10 (right) quickly converge as a function of N^. Using the convergence 



criteria described in ^64j , we choose A^^ = 3 as the optimal truncation order for the cumulant 
expansion in this example. 
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FIG. 9: Conventional {m^ff) and cumulant {rrv^r'l ) effective mass plots for A^ = 50 unitary 
fermions with Ar = 2 on a L = 12 lattice. Upper-left panel shows conventional, upper-right shows 
cumulant with N^ = 2, lower- left shows cumulant with A^^^ = 3 and lower-right shows cumulant 
with Nf^ = 4. The purple band in the effective mass plots represent fit results to the plateau region 
when one exists (in the case of Nf^ = 2, a constant plus an exponential fit function was used), the 
gray data in the cumulant effective mass plots represent the effective mass obtained by using the 
conventional method. 



To improve our estimates of the ground state energy, we have considered several different 



choices for the free parameter (3 appearing in the sink wavefunction defined in Eq. 13 The 



optimal values which we considered are provided in Table IIL By performing simultaneous 
fits to the correlators at different values of /3, we are able to obtain greater reliability in 
our ground state energies. As an example, we plot the effective masses for A^ = 60 unitary 
fermions on an L = 14 lattice using the cumulant expansion method truncated aX Nk = ^ 
for three values of (3 in Fig. Ill] (left). In the same figure (right) is a plot of the results from 
a simultaneous fit to all three correlation functions as a function of the starting time of the 
fitting region (the endpoint defined as rmax = 25 is held fixed). Using such simultaneous 



26 



0.45 
0.44 
gO.43 
H 0.42 
0.41 
0.40 



N=10 



Nk 




FIG. 10: Left: data points are the energies obtained using the cumulant effective mass for A^ = 10 
fermions, purple band is the energy obtained by using the conventional method. Right: data points 
are the energies obtained using the cumulant effective mass for A^ = 50 fermions. 




FIG. 11: (Left) Effective mass plot obtained with an N,^ = 4 truncation for A^ = 60 untrapped 
unitary fermions on a L = 14 lattice. Yellow diamonds, purple squares, and blue circles correspond 
to ensembles with sinks using /? = 1.0, 0.75, and 0.5, respectively. The dashed line represents 
the statistical uncertainty from a simultaneous fit, while the purple band represents the combined 
fitting statistical and systematic uncertainties. (Right) Results for simultaneous fits to the data in 
the left plot as a function of the beginning of the time interval used for fitting. The endpoint was 
held fixed at rmax = 25. 

fits, we are able to extract reliable energies from the plateau region where the three effective 
masses are statistically indistinguishable. 

The ground state energies for up to 66 unpolarized unitary fermions at finite volume were 
estimated using the analysis techniques described above and reported in Table [V] along with 
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TABLE V: Energies in units of the free gas energy for N = N^ + A^~'' paired unitary fermions 
in a finite box. Extrapolated results for L = oo reflect the fit parameter cq{N) obtained using a 
three-parameter fit described in the text. The uncertainties represent the fitting statistical and 
systematic uncertainties combined in a quadrature. 
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"^The ground state energies for TV = 4 are calculated from the ensembles used in Sec. IV C 



the truncation order A^^ used to obtain the result at each N. Energies are quoted in units of 
the thermodynamic hmit definition of the free gas energy, EpreeiN/V), for each value of A^ 
and L. The quoted errors represent both fitting statistical and systematic errors combined 
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FIG. 12: (Left) An infinite volume extrapolation of the ground state energy for A^ = 58 unitary 
fermions. (Right) Ground state energy extrapolated to zero density as a function of 1/A^. The red 
band represents a constant fit to the energies for 40 < A^ < 66. 



in quadrature as discussed in Sec. IV B For each fixed value of A^, we performed an infinite 



volume (or equivalently a zero density) extrapolation of the energy using data obtained at 
different volumes. As was the case for A^ = 4 fermions, we expect the leading and subleading 
volume dependence of the energies obtained using No = 5 tuned couplings to scale as L~^ 
and L~^-^^, respectively, corresponding to effects induced by £ = 1 two- and three-body 
operators. We therefore used the fit function Cq{N) + ci{N)/L^ + C2{N) / L^'^^ to perform an 
infinite volume the extrapolation. An example of such a fit for A^ = 58 untrapped fermions 



is shown in Fig. 12 (left). 



The infinite volume extrapolated energies E{N)/ EpreeiN) = cq{N) are tabulated in 
Table M and plotted in Fig. 12 (right) as a function of the inverse of fermion number. Our 
results show that the shell structure is present in the first and second shells (4 < A^ < 38), 
which is much more evident in the energies at finite volume. On the other hand, we find little 
evidence for shell effects within the last two shells (i.e., 40 < A^ < 66), suggesting that within 
the numerical uncertainty of our measurements, we are sufficiently near the thermodynamic 
limit to perform a thermodynamic limit extrapolation of the Bertsch parameter, given by 
^ = limAr_j.oo co(A^). Note that the A^-dependence is expected to be correlated since the 
energies at different A^ were determined from the same ensemble. To estimate the Bertsch 
parameter, we have performed a correlated constant fit to the infinite volume extrapolated 
energies over the fit range A^ G [40, 66], obtaining the estimated value: ^ = 0.3661q'q}^. 

The Bertsch parameter has been extensively studied in the past using quantum Monte- 
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Carlo (QMC) simulations. The earliest works based on a variational approach found an 
upper bound of .^ < 0.42(1) [2H1 EO], while a more recent QMC calculation for A^ = 66 with 
an extrapolation to zero range reported an upper bound of ^ < 0.383(3) [111112]. Numerous 
lattice simulations of two-component fermions in the unitary limit have been reported at both 
zero and nonzero temperature. References [39] and [38j quoted the Bertsch parameter values 
0.292(24) and 0.37(5), respectively, from finite temperature lattice simulations extrapolated 
to zero temperature. A different zero temperature lattice calculation with an infinite volume 
extrapolation for A^ = 10 and A^ = 14 yielded ^ = 0.292(12) and 0.329(5), respectively [37]. 
The Bertsch parameter has also been measured in several atomic experiments by studying 
pair correlation and absorption rates of ^Li and ^''K in a harmonic trap. Some recent 
experimental measurements reported ^ = 0.39(2) [lOj by Duke and 0.41(1) [11] by the Paris 
group. The most recent experimental determination by the group from the Massachusetts 



Institute of Technology (MIT) found 0.376(4) [13] . In Fig. 13, we summarize all analytical, 
numerical and experimental estimates of C, to date along with our value of the Bertsch 
parameter obtained from the simulations of up to A^ = 66 untrapped unitary fermions. 



References for the historical results are provided in Table VI Our determination of the 



Bertsch parameter appears as the latest data point in Fig. [13] and is statistically consistent 
with other recent findings. 

V. CONCLUSION 

We have studied up to 66 unpolarized unitary fermions in a periodic box by applying a 
lattice Monte Carlo method developed for studying large numbers of strongly interacting 
nonrelativistic spin-1/2 fermions ^50J. Our method differs from methods used in the past 
in that it does not make use of importance sampling, nor is it variational in nature. As 
such, our approach not only allows us to study unpolarized Fermi systems, but also systems 
with unequal numbers of spin up and spin down fermions. One of the main obstacles in 
calculating ground state energies of large numbers of fermions using our method is that 
it exhibits a severe distribution overlap problem, resulting in unreliable estimates of cor- 
relation functions. To solve this problem, we use a cumulant expansion technique for the 
logarithm of correlators [64] , which allows us to determine energies in a reliable manner with 
controlled systematic errors. The successful application of our method to unitary fermions 
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FIG. 13: Historical results for the Bertsch parameter determined experimentally, by analytic cal- 



culation, and by numerical simulation. Numerical values and citations are tabulated in Table VI 
our value is indicated as the latest simulation data point. 

gives us confidence that these techniques may prove useful in other situations where im- 
portance sampling is difficult. Conventional importance sampling schemes for Fermi gas 
calculations often use the A^-body correlator itself as an importance measure and so the en- 
semble generated is only of use for estimating a single observable for which it was designed. 
Our approach offers an advantage over such importance sampling schemes in that one may 
use the ensemble generated to reliably estimate all desired observables. Thus our approach 
avoids the multiplicative enhancement in computational cost by the number of measured 
observables which is inherent in calculations based on importance sampling. 
Our main findings for this study are summarized as follows: 

1. The exact diagonalization of the transfer matrix for two, three, and four fermions en- 
ables us to verify our simulation results, and to study systematic errors from spatial and 
temporal discretization. Our results for the spectrum of three and four fermions are 
in good agreement with benchmark calculations from other groups. While few-body 
systems were used in this paper as a way to test our methodology, they are interesting 
in their own right, and it looks feasible to use our methods in the future to measure 
the fascinating anomalous scaling behavior expected of three-body interactions. 
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TABLE VI: Historical results for the Bertsch parameter ^ determined experimentally (exp.), by 
numerical simulation (sim.) and by analytic calculation (anal.), along with publication (pub.) date. 
Values obtained variationally are upper bounds, and are indicated with an asterisk; simulation 
results without a quoted error bar should be regarded as approximate. 
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2. As part of our study of systematic effects, we fiave calculated the two-body s- and 
p-wave scattering phase phase shifts and few-body excited state energies in the lattice 
theory for various choices of L and No- 

3. Due to the highly improved, Galilean invariant action for which the first few terms in 
the effective range expansion for s-wave scattering have been systematically eliminated, 
we find mild volume dependence in the energies for the four volumes considered. The 
remaining finite- volume or discretization effects, where the leading contributions come 
from the p-wave and three-body operators, are eliminated by performing an infinite 
volume or equivalently a zero density extrapolation. 

4. The many-body ground state energies (measured in units of the noninteracting ener- 
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gies) show no discernible shell effects for A^ > 40, implying that the system is near the 
thermodynamic limit and therefore a reliable determination of the Bertsch parameter 
is possible. We determined the Bertsch parameter to approximately 4% statistical 
and systematic uncertainties and find agreement with the most recent experimental 
and numerical determinations by other research groups, thus demonstrating the suc- 
cess of our lattice construction and novel analysis methods as applied to many-body 
calculations. 

Our work shows that by combining novel statistical techniques and by perfecting the 
action it is possible to perform lattice Monte Carlo calculations that are quite competitive 
with methods that employ costly importance sampling. It is possible that the operator 
basis in Eq. |9]that we chose for perfecting the interactions is not the optimal set, and that 
our method for fixing the operator coefficients is not the optimal strategy. An interesting 
direction for future research would be to understand better whether there exists such an 
optimal strategy for perfecting the action, and whether there are benefits in combining both 
the perfect action and importance sampling techniques to further extend the computational 
reach of simulations for trapped atoms. 

Note added: After this work was completed, Ref. [7D] appeared using similar methods as 
described here, with compatible results for the unextrapolated Bertsch parameter at similar 
volumes. 
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Appendix A: Transfer matrices 

The ensemble average of a direct product of A^ propagators K~^(t;0) yields a V^- 
dimensional matrix: 

Un{t) = {K'\t; 0) ® . . . ® K-\t; 0)) , (Al) 

" V ' 

N 

which may be related to the Euclidean time-evolution operator for a system of A^ particles. 
Specifically one may define the time evolution operator as 

UN{r) = W,v(0)-^/'Wjv(r)W^(0)-i/2 , (A2) 

which satisfies the properties: 

Un{0) = 1 . (A3) 

One may then derive an analytic expression for the A^-particle transfer matrix, given by 
T/v = WAf(l), the eigenvalues of which yield the exponentiated energies of the iV-body 
system. 

Explicitly, the matrix elements of the A-body transfer matrix in momentum space is 
given by: 



W„...,<,'„|r„|<,,...,<,.>^n^^ 



N 



Il^<^ 



4 = 1 

N N 



(A4) 



where the ellipses represents higher order terms involving two or more contributions from 
the interaction C. In addition, the ellipses include contact terms which come from a slight 
modification of Wick's theorem in the case of Z2 fields. ^ Such contributions only appear in 
the case of four or more particles, however. 

In the case of A = A''^ + N-^ fermions, the above transfer matrix must be antisymmetrized 
with respect to momenta corresponding to each species. Although the transfer matrix is V^ 

^ For example, (^x'/'y'/'z'/'w) = (^xyiJzw + 'Jxz'^yw + '^xw'Jyz + w^xy'^yz'^zw, where w = for Gaussian and 
w = 1 for Z2 fields. 
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dimensional, we may reduce the dimensionality by a power of volume by projecting onto the 
center-of-mass (cm.) frame using the projection operator: 

N 

(q'l, • • • , q'jvl^c.m. |qi, . . . , q^v) = <5q;+...+q'^,o JJ ^<'^^ ' *^^^) 

A further reduction of the dimensionality may be performed by projecting the trans- 
fer matrix onto the positive and negative parity irreducible representations (irreps) r = 
Af , Af , E^ ,T^ ,T^ of the octahedral group Oh, using the projection operator: 

1 ^ 

(q'l, ... , q'^|Pr-|qi, • • • , qTv) = ^ 5^ Xr(^) n ^Ri9H,^^ ■ (A6) 

9 «=1 

In this expression, the sum is over the 48 group elements g of Oh, Xr{g) are the characters of 
the representation r, and R{g) are the three-dimensional rotation matrices corresponding to 
each group element g. The irreps r have dimensionality dr = 1, 1, 2, 3, and 3, respectively. 
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